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Abstract 

This  paper  describes  the  application  of  adaptive  multigrid  techniques  to  the  prob¬ 
lem  of  tropical  cyclone  track  prediction.  Based  on  the  nondivergent  barotropic  vorticity 
equation,  the  model  uses  an  adaptive  multigrid  method  to  refine  the  mesh  around  the 
moving  vortex.  Like  conventional  nested-grid  models,  this  model  achieves  nonuniform 
resolution  by  superimposing  uniform  grids  of  different  mesh  sizes.  Unlike  nested-grid 
models,  multigrid  processing  uses  the  interplay  between  solutions  on  fine  and  coarse 
grids — in  regions  where  they  overlap — to:  (1)  solve  the  implicit  problem  for  the  stream- 
function  with  optimum  efiiciency,  (2)  automatically  achieve  two-way  interaction  at  the 
grid  interfaces,  and  (3)  provide  accurate  truncation  error  estimates  for  use  in  deter¬ 
mining  where  to  refine  or  coarsen  the  grids.  An  exchange  rate  algorithm  accomplishes 
the  latter  task,  approximately  optimizing  the  grid  selection  based  on  a  user-specified 
tradeoff  between  accuracy  and  computational  work. 

Numerical  results  demonstrate  that  the  model  chooses  reasonable  grids  with  mini¬ 
mal  user  intervention.  Using  adaptive  mesh  refinement  is  at  least  an  order  of  magnitude 
more  efiicient  than  using  a  single  uniform  grid,  and  the  overhead  cost  of  adaptive  re- 
gridding  is  less  than  two  percent  of  the  total  execution  time.  The  adaptive  multigrid 
approach  allows  track  prediction  errors  due  to  discretization  to  be  essentially  eliminated 
from  the  problem  at  a  reasonable  computational  cost. 
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Introduction 


In  many  physical  problems  the  spatial  scale  of  the  solution  varies  significantly  over  the 
domain.  For  example,  a  tropical  cyclone  is  a  relatively  small-scale  vortex  embedded  in  a 
larger-scale  flow,  and  accurate  prediction  of  the  vortex  track  may  require  resolving  the  flow 
both  within  and  around  the  storm.  Solving  such  problems  numerically  often  requires  variable 
resolution  to  accurately  resolve  the  small-scale  features  without  compromising  the  overall 
efficiency.  Conventional  nested-grid  techniques  (Kurihara  et  al.,  1979;  Berger  and  Oliger, 
1984)  achieve  local  mesh  refinement  by  superimposing  nested  uniform  grids  of  varying  mesh 
sizes.  Such  methods  have  been  widely  used  for  many  years  (e.g.,  Kurihara  and  Bender, 
1980;  Juang  and  Hoke,  1992;  Skamarock  and  Klemp,  1993).  The  success  of  conventional 
nested-grid  tropical  cyclone  models  such  as  the  GFDL  (Kurihara  et  al.,  1998)  and  VICBAR 
(DeMaria  et  al.,  1992)  models  demonstrates  the  value  of  this  approach. 

In  contrast  to  nested-grid  models,  multigrid  methods  (Brandt,  1977)  use  multiple  grids 
of  different  mesh  size  covering  the  same  region.  The  goal  here  is  to  solve  the  discretized 
equations  with  optimum  efficiency:  using  simple  relaxation  on  each  grid  to  smooth  the  error 
on  the  scale  of  that  grid  produces  the  solution  (on  the  finest  grid)  to  the  level  of  truncation 
error  in  computational  work  proportional  to  the  number  of  unknowns  (usually  in  the  work 
of  just  a  few  relaxation  sweeps).  In  addition  to  speeding  up  the  solution  process,  multigrid 
processing  provides  accurate  truncation  error  estimates,  which  can  be  used  to  determine 
where  to  refine  (or  coarsen)  the  grids  in  a  local  refinement  scheme  or  to  achieve  higher-order 
accuracy  via  extrapolation. 

While  the  multigrid  literature  has  many  examples  of  combining  multigrid  processing  with 
nested  grids  (e.g.,  Ciesielski  et  al.,  1986;  Bai  and  Brandt,  1987;  Saleh,  1994;  Ruge  et  al., 
1995),  in  nested-grid  meteorological  models  the  interaction  between  the  solution  on  different 
grid  levels  has  in  general  not  been  exploited.  Indeed,  some  of  the  complexity  of  the  VICBAR 
model  is  due  to  the  fact  that  the  computational  grids  are  specifically  not  allowed  to  overlap. 

The  purpose  of  this  paper  is  to  demonstrate  the  potential  of  adaptive  multigrid  meth¬ 
ods  for  time-dependent  problems  requiring  local  mesh  refinement  by  describing  an  adaptive 
multigrid  tropical  cyclone  track  model.  In  the  tradition  of  the  barotropic  models  SANBAR 
(Sanders  et  al.,  1975)  and  VICBAR,  we  refer  to  this  multigrid  barotropic  model  as  MUD- 
BAR.  Our  primary  focus  here  is  on  the  numerical  method:  describing  the  adaptive  multigrid 
approach  in  detail  and  documenting  its  performance.  Indeed,  the  MUDBAR  model  currently 
has  the  simplest  possible  dynamics  (nondivergent  barotropic  equation)  and  no  physics.  Nev¬ 
ertheless,  the  model  is  fast  and  accurate,  and  could  be  useful  for  operational  forecasting  (in 
particular,  ensemble  forecasts)  if  augmented  with  boundary  data  from  a  global  model  and 
an  appropriate  initialization  scheme. 

The  remainder  of  the  paper  is  organized  as  follows.  Section  2  describes  the  model  equa¬ 
tions  and  solution  method.  The  grid  adaptation  algorithm  is  presented  in  section  3.  Numer¬ 
ical  results  appear  in  section  4,  and  our  conclusions  are  summarized  in  section  5. 
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2  Model  Description 


In  this  section  we  briefly  describe  the  governing  equations,  discretization,  grid  structure,  and 
solution  method  of  the  MUDBAR  model.  A  preliminary  version  of  the  model  was  described 
in  Fulton  (1997),  to  which  the  reader  is  referred  for  some  of  the  more  basic  details. 


2.1  Governing  Equations 


We  formulate  the  model  on  a  section  of  the  sphere,  transforming  longitude  A  and  latitude  (f> 
to  Cartesian  coordinates  x  and  y  via  the  Mercator  projection 

x  =  {\-  Xc)a  cos  (f>c,  y  =  [tanh~^(sin  <j))  -  tanh"^(sin  <l)c)]  a  cos  (pc,  (1) 

where  a  is  the  radius  of  the  earth,  so  the  projection  is  true  at  {KAc)  where  {x,y)  =  (0,0). 
The  model  consists  of  the  modified  barotropic  vorticity  equation 


29{i’X) 

d{x,y) 


+ 

ux 


(2) 


where  the  relative  vorticity  C  and  streamfunction  ip  are  related  by 

—  7^)  ^  =  C-  (^) 

Here  =  d'^/dx^  +  d^^/dy^,  /?  =  2Q,a~^  cos(p  (with  the  rotation  rate  of  the  earth),  and 
m  =  cos(pc/cos<p  is  the  map  factor.  The  /?-plane  approximation  is  recovered  by  setting 
m  =  1  and  P  =  2Qa~^  cos  (pc-  There  are  two  quasi-physical  parameters:  the  diffusion 
coefficient  u,  and  the  parameter  7  (inverse  of  the  effective  Rossby  radius)  which  helps  prevent 
retrogression  of  ultralong  Rossby  waves.  The  model  domain  is  a  rectangle  in  x  and  y  centered 
at  (x,y)  =  (0,0).  At  the  boundaries  we  specify  the  streamfunction  ip  (and  thus  the  normal 
component  of  the  velocity);  where  there  is  inflow,  we  also  specify  the  vorticity  With  these 
conditions  the  problem  is  well-posed  (Oliger  and  Sundstrom,  1978). 


2.2  Discretization 

The  space  discretization  uses  second-order  finite  differences  on  uniform  rectangular  grids  (as 
detailed  below),  approximating  the  advection  terms  by  the  Arakawa  Jacobian  (Arakawa, 
1966).  Details  of  the  boundary  discretization  are  given  in  Fulton  (1997).  On  a  uniform  grid 
with  mesh  size  h  in  x  and  y,  the  space-discretized  approximations  to  (2)  and  (3)  can  be 
represented  in  the  form 

^  =  F'‘(V’",C'*)  (4) 

and 

Ihiph  ^  (5) 
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where  the  grid  functions  and  C'*  consist  of  the  values  of  the  approximate  solution  on  the 
grid  h. 

The  time  discretization  uses  the  classical  fourth-order  Runge-Kutta  (RK4)  scheme;  this 
is  highly  accurate,  allows  relatively  large  time  steps  for  stability,  and — since  it  is  a  one- 
step  scheme — is  easy  to  implement  and  has  no  spurious  computational  modes.  This  scheme 
involves  four  stages  per  time  step,  each  of  which  predicts  a  new  value  for  C'*  via  (4)  and 
solves  (5)  for  the  corresponding 


2.3  Grid  Structure  and  Local  Time  Stepping 

The  above  discretizations  are  embedded  in  an  adaptive  method  which  includes  local  refine¬ 
ments  in  both  space  and  time,  superimposing  nested  uniform  grids  with  different  mesh  sizes 
to  adapt  the  resolution  near  the  storm.  The  base  grid  Gi  with  mesh  size  hi  covers  the  entire 
computational  domain,  while  successively  finer  patches  Gi  with  mesh  sizes  hi  {I  =  2,  , 

n)  cover  smaller  nested  areas  as  shown  in  the  example  in  Fig.  1.  We  use  the  mesh  ratio 
hi-i/hi  =  2  (to  facilitate  multigrid  processing),  and  require  the  patches  to  be  strictly  nested 
(fine  grid  contained  in  the  interior  of  the  next  coarser  grid)  and  aligned  (patch  boundaries 
coincide  with  coarse-grid  lines).  For  simplicity,  we  allow  only  one  region  of  refinement  (i.e., 
one  grid  patch  per  mesh  size),  but  this  is  not  a  fundamental  limitation  of  the  method. 

The  time  stepping  algorithm  uses  local  time  stepping,  i.e.,  smaller  time  steps  on  finer 
grids;  this  algorithm  is  similar  to  that  used  in  most  nested-grid  models  (e.g.,  Berger  and 
Oliger,  1984;  DeMaria  et  ah,  1992).  Specifically,  with  two  grids  (coarse  and  fine)  one  full 
time  step  is  executed  as  follows: 

1.  One  step  (length  At)  on  the  coarse  grid, 

2.  Two  steps  (length  At/2)  on  the  fine  grid,  using  boundary  values  interpolated  from 
the  coarse  grid  in  space  (C'*  linearly  and  by  cubic  interpolation)  and  time  (linear 
interpolation), 

3.  Transfer  the  fine-grid  solution  to  the  coarse  grid  where  they  overlap  by  injection 
and  by  full  weighting). 

This  algorithm  generalizes  recursively  to  more  than  two  grids.  Its  application  to  three 
computational  grids  is  illustrated  in  Fig.  2. 

2.4  Multigrid  solution  for  streamfunction 

To  solve  (5)  eflficiently  for  the  streamfunction  on  any  computational  grid  Gi  we  use  a 
multigrid  method.  This  uses  point  Gauss-Seidel  relaxation  with  red-black  ordering,  full 
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Figure  1;  Sample  computational  grids  (h  =  128,  64,  32,  16km). 


weighting  of  residuals,  bilinear  interpolation  of  corrections,  and  the  full  multigrid  (FMG) 
algorithm  with  one  V(l,l)-cycle  per  level  and  bicubic  initial  interpolation.  As  these  elements 
are  all  standard,  the  reader  is  referred  to  Fulton  et  al.  (1986)  and  the  references  therein  for 
details.  The  resulting  algorithm  solves  accurately  for  the  streamfunction,  producing  residuals 
smaller  than  the  truncation  error;  the  computational  work  is  comparable  to  about  eight 
relaxation  sweeps  on  the  finest  grid,  and  overall  accounts  for  about  two-thirds  of  the  total 
execution  time  of  the  model^.  Numerical  results  show  that  using  more  sweeps  per  cycle  (or 
more  cycles)  produces  smaller  residuals  but  no  significant  improvement  in  the  track  error, 
and  thus  is  not  worth  the  extra  computational  expense. 

When  used  in  the  context  of  the  adaptive  grid  structure  described  above,  the  standard 
multigrid  approach  is  modified  in  several  ways  which  deserve  mention.  First,  with  local  time 
stepping,  the  next  coarser  grid  Gj_i  may  be  at  a  different  time  t  (e.g.,  at  the  end  of  time 
steps  2,  3,  and  6  in  Fig.  2).  In  this  situation  the  main  computational  grids  cannot  be  used 
for  multigrid  processing.  Instead,  for  each  computational  grid  Gi  we  use  an  additional  set 

^Achieving  the  same  accuracy  with  single-grid  relaxation  (SOR  with  optimal  relaxation  parameter)  re¬ 
quires  about  50-100  sweeps  and  roughly  ten  times  as  much  total  computer  time. 
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Figure  2:  Example  of  time  stepping  algorithm  for  one  full  time  step  At  on 
three  computational  grids.  The  steps  are  executed  in  the  numbered  order, 
with  left  arrows  denoting  the  transfer  of  the  solution  from  the  fine  grid  to  the 
coarse  grid  where  they  overlap.  Solid  circles  mark  times  and  grids  where  the 
truncation  error  is  estimated  for  regridding. 


of  coarse  grids  with  mesh  sizes  2hi,  4hi,  Shi,  etc.,  covering  the  same  domain  as  Gi  (Pulton, 
1997).  These  “local  coarse  grids”  may  be  used  in  two  ways.  The  Berger-Oliger  (BO)  method 
uses  them  at  all  times,  resulting  in  only  one-way  interaction  between  the  grids:  the  fine  grid 
uses  boundary  values  from  the  coarse  grid  (step  2  above),  but  its  effect  is  felt  on  the  coarse 
grid  only  through  the  solution  transfer  in  the  overlap  region  at  the  end  of  the  time  step 
(step  3  above).  In  contrast,  the  multigrid  (MG)  method  uses  the  full  computational  grid 
Gi-i  instead  of  the  first  local  coarse  grid  during  multigrid  cycling  in  the  final  RK4  stage  of 
the  second  fine-grid  time  step  (e.g.,  at  the  end  of  time  steps  4,  5,  and  7  in  Fig.  2).  This 
results  in  automatic  two-way  interaction  between  the  computational  grids,  since  fine-grid 
information  can  affect  the  entire  coarse  grid  immediately  through  the  relaxation  process. 
Preliminary  results  (Fulton,  1997)  showed  that  the  differences  between  these  two  methods 
are  generally  small;  further  experimentation  (see  section  4.4)  shows  some  improvement  with 
the  MG  approach,  so  that  is  used  in  all  other  results  presented  here. 

Second,  using  the  full  computational  coarse  grid  in  the  MG  method  requires  Pull  Ap¬ 
proximation  Scheme  (FAS)  processing  (Brandt,  1977),  since  each  grid  patch  covers  a  smaller 
domain  than  the  previous  one.  Specifically,  if  on  a  grid  patch  Gi  (/  >  1)  with  mesh  size 
h  =  hi  -we  denote  by  the  current  approximation  to  the  true  (discrete)  solution  ip^,  then 
the  diagnostic  equation  to  be  solved  in  the  overlap  region  on  the  next  coarser  grid  Gi-i  with 
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mesh  size  2h  =  is 


l2h^2h  ^  ^2h  ^  j2h(^^h  _  (5) 

and  the  corresponding  correction  to  the  fine-grid  solution  is  given  by 

+  (7) 

Here,  the  fine-to-coarse  transfer  operators  and  represent  injection  and  full  weighting, 
respectively,  and  the  coarse-to-fine  transfer  operator  represents  bilinear  interpolation.  In 
the  region  of  the  coarse  grid  not  covered  by  the  fine  grid,  the  FAS  equation  (6)  is  replaced 
by  and  no  correction  (7)  is  needed.  When  the  local  coarse  grids  are  used  (e.g., 

for  the  first  time  step  on  the  fine  grid)  FAS  processing  is  still  used  for  convenience,  even 
though  the  simpler  correction  scheme  could  be  used. 


Third,  to  properly  represent  the  net  fine-grid  vorticity  on  the  coarse  grid  when  the  grids 
are  nested,  a  correction  must  be  applied  at  the  grid  interfaces.  As  shown  by  Bai  and  Brandt 
(1987),  this  correction  is  accomplished  by  setting  the  coarse-grid  vorticity  at  the  grid  interface 

(S) 


2h 


2h 


Here  Bl^  is  the  one-dimensional  full  weighting  operator  along  the  boundary  and  is  the 
one-sided  difference  approximation  (over  grid  length  2h)  to  the  outward  normal  derivative^. 
For  example,  at  a  point  along  the  eastern  boundary  of  a  fine-grid  patch  indexed  in  x  and  y 
as  (6,y)  on  the  fine  grid  and  (B,  J)  on  the  coarse  grid,  we  have 


{Dn  i’  hj  -  - 2h - ’  ^  "  2h 


(9) 


and 


(BfC'*) 


— 

B,J 


J_  J_ 


(10) 


where  is  the  initial  FAS  approximation  on  the  coarse  grid.  Other  than  the  cor¬ 

rection  (8),  no  special  treatment  is  required  at  the  grid  interfaces;  in  particular,  no  extended 
overlap  or  transition  region  is  used  or  needed. 

^In  Bai  and  Brandt  (1987)  the  sign  of  D^'*  is  reversed,  which  appears  to  be  an  error. 
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The  Grid  Adaptation  Algorithm 


The  preliminary  version  of  the  model  described  in  Fulton  (1997)  used  movable  grids,  i.e., 
patches  with  sizes  fixed  in  advance  were  simply  moved  when  necessary  to  keep  them  ap¬ 
proximately  centered  on  the  vortex.  The  results  indicated  that  some  choices  of  patch  sizes 
substantially  increased  the  accuracy  or  efficiency  (relative  to  using  a  single  uniform  grid), 
while  others  led  to  significantly  less  improvement.  Since  the  optimal  combination  of  patch 
sizes  will  depend  on  both  the  vortex  and  the  surrounding  flow  and  may  vary  as  the  solution 
evolves,  choosing  fixed  patch  sizes  “manually”  is  less  than  ideal:  a  fully  automatic  algorithm 
would  be  advantageous.  This  section  details  such  an  algorithm;  we  refer  to  the  grids  it 
chooses  as  adaptive  grids. 


3.1  Basic  regrid  strategy 

The  MUDBAR  model  regrids,  i.e.,  chooses  patch  sizes  and  locations  adaptively,  by  estimating 
the  truncation  error  and  placing  the  grid  patches  where  the  truncation  error  is  large.  The 
regridding  algorithm  is  embedded  in  the  local  time  stepping  algorithm  of  section  2.3  as  follows 
(see  Fig.  3).  Before  starting  a  time  step  on  a  given  grid  with  mesh  size  h — here  considered  to 
be  the  “coarse  grid” — the  truncation  error  is  estimated  (see  section  3.2).  Points  where  the 
truncation  error  is  large  (see  section  3.3)  are  flagged  for  refinement,  as  indicated  by  “x”  in 
Fig.  3(a).  Second,  after  the  the  time  step  is  complete  the  truncation  error  is  again  estimated 
and  large  values  are  flagged  for  refinement,  as  shown  in  Fig.  3(b).  Third,  the  smallest 
acceptable  patches  at  the  beginning  and  end  of  the  step  are  determined  (dotted  and  dashed 
outlines,  respectively,  in  Fig.  3),  and  a  suitable  fine-grid  patch  with  mesh  size  h/2  is  chosen 
to  cover  both  of  these  (solid  outlines  in  Fig.  3(b)).  Finally,  that  patch  is  “filled”  with  values 
of  the  vorticity  and  streamfunction  The  patch  so  constructed  is  then  used  for 

the  two  fine-grid  time  steps  which  coincide  with  the  single  time  step  originally  taken  on  the 
coarse  grid. 

This  process  is  repeated  recursively,  i.e.,  during  each  of  the  two  fine-grid  steps  a  size 
and  location  are  chosen  for  the  next  finer  patch  (if  needed),  proceeding  to  increasingly  finer 
grids  until  either  further  refinement  is  not  needed  or  a  specified  maximum  number  of  grids 
is  reached.  This  recursion  is  illustrated  in  Fig.  2,  where  the  solid  circles  denote  times  at 
which  the  truncation  error  is  estimated.  For  example,  the  truncation  error  on  level  /  =  2  is 
computed  before  and  after  time  step  5  (at  t  =  At/2  and  t  =  At,  respectively)  and  used  to 
determine  the  location  and  size  of  the  next  finer  patch  Z  =  3  for  time  steps  6  and  7  (from 
t  =  At/2  to  t  =  At). 

To  fill  a  given  patch,  the  vorticity  is  set  by  copying  values  from  the  previous  version 
of  that  grid  (if  any)  where  they  coincide  and  using  bilinear  interpolation  from  the  coarse 
grid  where  the  fine  grid  is  new;  at  t  =  0  the  specified  initial  vorticity  is  used  instead.  The 
streamfunction  is  set  likewise  except  that  bicubic  interpolation  is  used;  if  the  patch 
is  new  (e.g.,  at  t  =  0)  then  is  obtained  instead  by  solving  (5)  with  boundary  values 
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Figure  3:  Choice  of  a  fine-grid  patch.  The  contours  show  truncation  errors 
(large  values  flagged  by  x)  at:  (a)  the  start  of  the  time  step,  and  (b)  the  end 
of  the  time  step.  The  boundaries  of  the  minimum  acceptable  patch  at  the 
start  (dotted)  and  end  (dashed)  are  shown,  along  with  the  patch  chosen  for 
the  step  (solid). 


interpolated  (cubically)  from  the  coarse  grid.  During  model  initialization  at  t  =  0,  adaptive 
grids  are  constructed  by  starting  with  a  specified  base  grid  and  constructing  patches  as 
dictated  by  the  truncation  error  as  outlined  above.  If  the  first  patch  covers  most  of  the  base 
grid  (i.e.,  the  base  grid  resolution  is  inadequate),  this  patch  is  extended  to  the  full  domain  to 
replace  the  base  grid  and  the  process  is  repeated.  Once  the  initial  grids  are  set — and  thus  the 
fine-grid  information  has  influenced  the  solution  on  all  coarser  levels-a  second  pass  through 
the  grid  adaptation  algorithm  is  made  to  take  advantage  of  that  fine-grid  information. 

Even  though  new  grids  are  chosen  for  each  time  step  on  each  level,  the  regridding  process 
adds  surprisingly  little  overhead  to  the  execution  time.  All  calculations  required  for  regrid¬ 
ding  (estimating  the  truncation  error,  choosing  the  grid,  and  filling  it  by  copying  and/or 
interpolation)  typically  amount  to  less  than  two  percent  of  the  total  execution  time  of  the 
model.  This  slight  overhead  is  more  than  made  up  for  by  the  fact  that  the  grid  chosen  at 
each  step  on  each  level  is  precisely  the  size  needed  for  the  entire  step,  and  no  larger.  Indeed, 
a  simpler  scheme  based  on  regridding  on  all  levels  only  at  the  end  of  each  base-grid  time 
step  regrids  much  less  often  but  must  choose  larger  patches  (since  they  must  be  adequate 
for  several  fine-grid  time  steps)  and  thus  takes  more  time  (typically  20-30%  more). 

The  model  also  includes  an  option  for  using  movable  grids,  i.e.,  fixed  patch  sizes  specified 
in  advance.  This  differs  from  the  preliminary  version  described  in  Pulton  (1997)  in  that 
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the  patch  locations  are  now  chosen  for  an  interval  in  time,  rather  than  an  instant.  This  is 
done  by  finding  the  vortex  center  at  the  beginning  and  end  of  the  coarse-grid  time  step  and 
centering  the  fine-grid  patch  for  the  corresponding  two  fine-grid  time  steps  halfway  between 
these  two  points.  Except  for  this  choice  of  patch  location  (and  the  fixed  size),  the  regrid 
strategy  with  movable  grids  is  identical  to  that  described  above  for  adaptive  grids. 


3.2  Truncation  error  estimates 


The  truncation  error  estimate  needed  for  grid  adaptation  may  be  computed  essentially  for 
free  using  FAS  processing;  this  idea  was  introduced  by  Brandt  (1977)  and  is  similar  to  that 
described  by  Skamarock  (1989).  Specifically,  we  want  to  base  our  grid  refinement  criterion 
on  the  truncation  error 

rh 

where  ip  and  C  represent  the  true  solution  of  the  continuous  problem  (3)  and  represents 
their  pointwise  restriction  to  grid  h.  It  can  be  shown  that  the  relative  truncation  error 

rf  :=  C'*  (12) 


provides  an  accurate  approximation  to  the  truncation  error  difference  In  contrast 

to  the  analysis  of  Bernert  (1997),  we  find  that  with  injection  as  the  fine-to-coarse  transfer 
in  (12),  this  estimate  is  fourth-order  accurate  (e.g.,  Fulton,  1989).  Since  r*  =  0{h‘^),  we 
obtain  the  estimate  ^ 

(13) 


,h  ^  t  2h 


which  has  accuracy  0{h^).  Since  is  known  only  at  the  coarse-grid  points,  it  is  possible 
to  compute  only  at  the  coarse-grid  points;  however,  the  estimate  (13)  is  scaled  so  it 
is  indeed  an  estimate  of  the  truncation  error  on  the  fine  grid.  Note  that  while  can  be 
computed  from  (12)  with  essentially  no  work  (only  one  subtraction  per  coarse-grid  point) 
during  the  last  V-cycle  of  the  FMG  method,  computing  it  separately  after  completing  that 
cycle  adds  very  little  work  and  produces  a  better  estimate  of  r^,  since  it  is  based  on  the  final 
solution  for  ip^  (rather  than  the  initial  approximation) . 


3.3  Patch  location  and  size 

The  determination  of  what  constitutes  a  “large”  value  of  truncation  error  follows  the  method 
outlined  by  Brandt  (1977),  which  uses  a  parameter  A,  called  the  exchange  rate,  to  control  the 
trade-off  between  increased  accuracy  and  increased  work®.  A  point  on  the  grid  h  is  flagged 
for  refinement  if 

hV'*  >  A  (14) 

^The  exchange  rate  A  originates  as  a  Lagrange  multiplier  in  the  problem  of  minimizing  the  error  subject 
to  the  constraint  of  constant  work  (or  minimizing  the  work  subject  to  constant  error). 
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there.  In  this  criterion,  accounts  for  the  work  needed  on  the  new  grid  and  r'*  measures  the 
accuracy.  Thus,  using  a  large  value  of  A  says  to  refine  only  when  the  error  is  large,  and  thus 
should  give  a  fast  solution — with  a  relatively  large  error.  Conversely,  using  a  small  value  of 
A  indicates  a  willingness  to  exchange  more  work  to  get  more  accuracy  (smaller  error). 

Once  the  points  where  the  truncation  error  on  the  coarse  grid  (mesh  spacing  h)  is  large — 
in  the  sense  of  (14)— are  located,  the  boundaries  of  the  fine-grid  patch  (mesh  spacing  h/2) 
must  be  determined.  For  simplicity,  the  model  currently  uses  only  rectangular  grids  (and 
only  one  grid  per  value  of  h,  so  only  one  region  is  refined),  so  we  can  simply  determine 
the  smallest  rectangles  enclosing  the  large  truncation  errors  before  and  after  the  step,  as 
shown  in  Fig.  3.  The  tentative  boundary  of  the  patch  for  the  full  time  step  is  simply  the 
smallest  rectangle  enclosing  both  of  these  (i.e.,  their  convex  hull).  This  boundary  may  then 
be  adjusted  by  adding  a  “buffer”  of  several  grid  intervals  (usually  two)  to  give  some  flexibility 
when  the  truncation  error  is  highly  localized.  Finally,  the  patch  size  may  be  increased  (or 
possibly  slightly  decreased)  if  needed  to  achieve  good  coarsenability,  i.e.,  so  that  the  number 
of  points  on  the  coarsest  local  coarse  grid  is  small,  thus  increasing  the  efficiency  of  solving  for 
the  streamfunction  on  the  patch.  The  algorithm  employed  here  uses  simple  work  estimates 
as  detailed  in  the  Appendix.  In  any  case,  strict  nesting  of  the  patch  is  enforced. 
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4  Results 


This  section  summarizes  the  results  of  many  model  runs  designed  to  explore  the  performance 
of  the  model  for  various  situations. 


4.1  Initial  conditions 


Except  as  otherwise  noted,  we  use  the  initial  conditions  of  DeMaria  (1985)  and  associated 
parameter  values  as  follows.  The  initial  vortex  has  tangential  wind  given  by 


KM  =  2K.(^) 


exp[-a(r/r^)^] 


(15) 


where  r  =  [{x  -  a^o)^  +  (v  -  is  the  radial  distance  from  the  vortex  center  {xo,yo)- 

Note  that  V  has  the  approximate  maximum  value  near  r  =  Vm  (exact  when  a  =  0); 
the  exponential  factor  is  included  to  make  V  vanish  quickly  for  large  r.  The  vortex  can 
be  made  slightly  elliptical  (to  introduce  some  small-scale  structure,  as  in  the  spiral  bands 
studied  by  Guinn  and  Schubert  (1993))  by  dividing  the  x  or  y  factor  in  the  definition  of  r 
by  (1  -  e^),  where  e  is  the  desired  ellipticity.  We  consider  two  cases:  a  “weak  hurricane” 
with  Vm  =  30ms“^  and  Vm  =  80  km,  and  a  “strong  hurricane”  with  Vm  =  60  ms"^  and 
Tm  =  20  km.  For  both  we  use  the  values  a  =  10"®,  6  =  6,  and  e  =  0.8.  The  environmental 
flow  is  the  zonal  current  with  streamfunction  given  by 

=  (It)  Yt)  ■ 

The  computational  domain  is  a  square  of  side  length  4096  km  on  a  /?-plane  centered  at 
latitude  20°  N.  The  vortex  is  centered  initially  at  xq  =  768  km  and  yo  =  —768  km,  and  we 
use  the  values  So  =  10  ms“\  L  =  4096  km,  and  7  =  0  and  i/  =  0. 


4.2  Sample  solution 

Figure  4  illustrates  the  solution  of  the  model  for  the  weak  hurricane  case.  The  left-hand  pan¬ 
els  show  the  streamfunction  ip  and  the  right-hand  panels  show  the  corresponding  vorticity  C- 
For  this  run,  the  base  grid  has  size  64  x  64  [hi  =  64  km),  the  exchange  rate  is  A  =  1000,  and 
the  grid  refinement  is  limited  to  three  patches,  so  the  finest  mesh  used  is  ^4  =  8  km.  The 
boundaries  of  the  grid  patches  chosen  by  the  adaptive  algorithm  are  shown  in  the  figures.  It 
can  be  seen  that  the  patch  sizes  tend  to  grow  slightly  during  the  model  run;  the  apparent 
slight  bias  toward  refining  the  grids  behind  the  vortex  is  likely  due  to  the  larger  truncation 
errors  in  the  wake  which  are  characteristic  of  centered  finite  differences.  This  eflfect  could 
probably  be  lessened  by  including  a  small  amount  of  dissipation;  however,  the  model  run 
shows  that  no  dissipation  is  needed  to  achieve  reasonable  performance.  In  particular,  no 
problems  are  observed  at  the  grid  interfaces. 
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4.3  Accuracy  and  efficiency 


To  measure  the  accuracy  of  the  model  in  the  absence  of  an  analytical  solution,  we  compare 
the  tropical  cyclone  track  to  that  obtained  in  a  high-resolution  reference  run  with  uniform 
resolution  h  =  4  km.  The  track  error  is  defined  as  the  distance  between  the  location  of  the 
vortex  center  (vorticity  maximum)  and  that  of  the  reference  run,  with  the  center  located  by 
quadratic  interpolation  using  four  points  surrounding  the  vorticity  maximum  on  the  finest 
grid  (Smith  et  ah,  1990);  the  mean  forecast  error  is  computed  over  a  72  hour  model  run  (using 
trapezoidal  quadrature  with  a  step  size  of  one  hour).  Computational  work  is  measured  by 
the  CPU  time  required  for  the  model  run  (on  a  SUN  UltrafiO  workstation  with  one  360  MHz 
processor). 

Figure  5  shows  the  mean  forecast  error  as  a  function  of  the  computer  time  used  for  a 
large  number  of  model  runs  for  the  weak  hurricane  case.  The  points  plotted  with  “-I-”  are  for 
uniform-resolution  model  runs  (i.e.,  a  base  grid  with  no  finer  grid  patches)  with  the  indicated 
resolution.  The  line  joining  these  points  indicates  a  baseline  (uniform  grid)  performance.  The 
isolated  points  to  the  left  of  this  line  are  for  model  runs  with  various  combinations  of  movable 
grids  (i.e.,  fixed  size  patches),  using  essentially  every  reasonable  combination  of  patches  with 
mesh  sizes  from  h  =  128km  ioh  =  4km.  The  points  plotted  with  “x”  are  for  the  adaptive 
model.  These  are  connected  by  lines  according  to  the  finest  mesh  size  allowed  in  the  run; 
different  points  on  each  line  correspond  to  different  values  of  the  exchange  rate  A  from  10^  to 
10^,  with  errors  generally  decreasing — and  CPU  time  increasing — as  A  decreases.  The  figure 
shows  a  savings  of  roughly  a  factor  of  ten  in  computer  time  by  using  local  mesh  refinement, 
and  suggests  that  in  most  cases  the  grids  chosen  by  the  adaptive  algorithm  are  reasonable. 


Figure  6  shows  analogous  results  for  the  strong  hurricane  case.  Since  the  vortex  is  both 
smaller  and  stronger,  errors  for  a  given  amount  of  work  are  generally  larger  than  in  the 
previous  case.  Again  the  adaptive  multigrid  method  chooses  reasonable  grids;  in  some  cases 
the  savings  of  execution  time  (compared  to  using  uniform  resolution)  approach  a  factor 
of  100.  Similar  results  are  obtained  in  other  cases.  For  example,  for  Figure  7  the  strong 
hurricane  vortex  is  embedded  in  a  different  environmental  flow  with  zonal  wind 


u{y)  = 


exp 


(17) 


where  uq  =  10ms“^  and  yi  =  443.4  km.  For  this  case  the  vortex  is  centered  initially  at 
xo  =  768  km  and  yo  =  768  km,  where  the  positive  Laplacian  of  the  environmental  vorticity 
leads  to  increased  sensitivity  of  the  vortex  track  to  initial  position  errors  (DeMaria,  1985). 
For  the  small  vortex  used  here,  this  sensitivity  is  not  evident;  the  adaptive  model  performs 
as  well  or  better  than  in  the  other  cases. 


4.4  Effect  of  two-way  interaction 

Figure  8  demonstrates  the  effect  of  including  two-way  interaction  between  the  computational 
grids.  Each  point  represents  the  difference  in  mean  forecast  error  between  the  Berger-Oliger 
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Figure  7:  Same  as  Fig.  6  except  for  the  environmental  flow  given  by  (17) 
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MEAN  ERRORS 


Figure  8:  Comparison  of  accuracy  achieved  by  the  BO  (one-way  interaction) 
and  MG  (two-way  interaction)  methods  for  the  weak  hurricane  case. 


(BO)  method  and  multigrid  (MG)  method  as  described  in  section  2.4.  The  points  are  plotted 
as  a  function  of  the  error  produced  by  the  MG  method  for  many  runs  with  different  choices 
of  movable  grids.  The  differences  between  the  accuracy  of  the  two  methods  are  small — and 
can  go  either  way — but  when  one  method  is  better,  it  is  usually  the  multigrid  method,  which 
includes  two-way  interaction. 


4.5  A  barotropically  unstable  vortex 

The  adaptive  multigrid  method  might  be  especially  advantageous  for  modeling  a  vortex 
which  not  only  moves  but  has  small-scale  structure  which  evolves  with  time.  While  the  non- 
divergent  barotropic  model  employed  here  is  incapable  of  representing  the  baroclinic  effects 
and  convection  related  to  tropical  cyclone  intensification,  it  can  represent  the  basic  physical 
process  of  potential  vorticity  (PV)  mixing.  Thus,  we  consider  the  case  of  a  barotropically 
unstable  vortex  which  develops  small-scale  asymmetries  due  to  chaotic  nonlinear  PV  mixing 
and  evolves  toward  a  smooth  symmetric  vortex  as  it  moves. 
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As  an  initial  condition  we  use  the  “PV  ring”  vortex  of  Schubert  et  al.  (1999)  embedded 
in  the  zonal  current  (16).  All  parameters  have  the  values  given  in  that  paper,  except  for 
the  diffusion  parameter  v  for  which  we  use  the  larger  value  300m^s"\  corresponding  to  an 
e-folding  time  of  338  s  for  a  wavelength  of  2  km.  For  this  case  we  run  the  model  with  spherical 
geometry  using  the  Mercator  projection  (1).  We  present  results  of  two  model  runs:  one  using 
movable  grids  and  the  other  using  adaptive  grid  patches  (with  exchange  rate  A  =  20).  In 
both  cases  we  start  with  the  base  grid  size  128  x  128  with  mesh  spacing  h  =  32  km  and 
corresponding  time  step  At  =  20  minutes.  The  movable  grid  run  uses  six  square  patches 
with  side  lengths  ranging  from  1/4  to  1/32  of  the  domain  size.  The  adaptive  run  is  limited 
to  six  adaptive  grid  patches;  however,  the  model  rejects  the  base  grid  as  inadequate  and 
consequently  uses  the  base  grid  mesh  spacing  h  =  16  km  and  up  to  five  adaptive  grids. 
Thus,  for  both  runs  the  model  domain  is  4096  km  square  and  the  finest  mesh  spacing  is 
h  =  0.5  km. 

Figure  9  shows  the  details  of  the  solution  computed  using  movable  grids  (left  panels)  and 
adaptive  grids  (right  panels).  Only  a  region  256  km  square  centered  on  the  vortex  is  shown. 
Since  the  vorticity  is  not  always  maximum  at  the  center  of  the  ring,  we  define  the  vortex 
center  as  the  vorticity  centroid  (taken  over  the  region  where  C  >  /,  to  eliminate  bias  from  the 
environmental  vorticity).  These  results  are  quite  similar  to  those  of  Schubert  et  al.  (1999)  up 
to  at  least  t  =  6  hours,  verifying  that  the  characteristic  wavenumber-four  pattern  observed 
in  that  study  is  due  to  the  dynamics,  and  is  not  an  artifact  of  the  periodicity  required  for 
the  spectral  model  used.  At  later  times,  the  solutions  are  qualitatively  similar  but  differ  in 
detail,  as  expected  due  to  the  chaotic  nature  of  the  PV  mixing.  Comparing  the  two  solutions, 
we  find  that  the  details  of  the  small-scale  structure  are  quite  similar;  also,  the  vortex  tracks 
(not  shown)  differ  by  only  about  7  km  over  the  course  of  the  72  hour  model  run.  However, 
the  run  using  adaptive  grids  was  six  times  faster  (about  16  minutes  on  a  SUN  Ultra60),  since 
as  the  solution  evolves  toward  a  smooth  symmetric  vortex  the  finest  grid  patches  are  no 
longer  needed  and  are  automatically  discarded. 
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Figure  9:  Vorticity  for  the  PV  ring  case  computed  using  movable  grids  (left 
panels)  and  adaptive  grids  (right  panels).  Solid  rectangles  show  patch  bound¬ 
aries.  The  finest  patch  shown  has  mesh  size  h  =  0.5  km  in  each  case  except 
the  adaptive  grids  at  t  =  24  hours  (h  =  1  km)  and  <  =  72  hours  {h  =  2  km). 
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5  Conclusions  and  Remarks 


We  have  used  the  problem  of  tropical  cyclone  track  prediction  to  investigate  the  potential  of 
the  adaptive  multigrid  method.  Combining  nesting  of  uniform  grids,  multigrid  processing, 
and  adaptive  mesh  refinement  based  on  truncation  error  estimates,  the  method  provides  a 
seamless  discretization  with  enhanced  resolution  precisely  where  it  is  needed.  Numerical 
results  from  the  nondivergent  barotropic  model  MUDBAR  show  the  method  performs  well, 
generally  reducing  the  execution  time  by  an  order  of  magnitude  or  more  (relative  to  using 
uniform  grids)  and  choosing  reasonable  grids  with  a  minimum  of  user  intervention.  Compared 
to  conventional  nested-grid  models,  the  treatment  of  grid  interfaces  is  particularly  simple 
(although  that  simplicity  may  be  more  a  result  of  using  simple  dynamics  rather  than  the 
numerical  method).  These  results  are  encouraging,  and  suggest  that  adaptive  multigrid 
methods  may  be  useful  in  other  problems  requiring  time-dependent  local  mesh  refinement. 

While  the  model  as  it  stands  is  fast  and  accurate,  there  is  much  room  for  improve¬ 
ment.  We  are  currently  working  on  two  modifications  which  should  substantially  improve 
the  inherent  accuracy,  namely,  converting  the  model  to  the  shallow-water  equations  and 
implementing  fourth-order  space  differencing.  To  be  useful  for  actual  prediction,  the  model 
should  be  supplied  with  time-dependent  boundary  data  from  a  global  model  and  an  appro¬ 
priate  initialization  scheme.  Given  its  speed  and  accuracy,  the  resulting  model  could  be 
especially  appropriate  for  ensemble  forecasting  techniques. 

The  numerical  method  presented  here  may  in  fact  be  overkill:  for  tropical  cyclone  track 
prediction,  discretization  errors  may  be  less  important  than  observational  uncertainties  or 
incomplete  representations  of  the  physics.  Indeed,  since  the  discretization  is  the  only  source 
of  error  considered  here,  our  numerical  results  show  far  more  accuracy  than  could  be  hoped 
for  in  operational  prediction.  Nevertheless,  the  adaptive  multigrid  approach  allows  one  to 
effectively  remove  discretization  error  from  the  problem,  and  does  so  in  a  natural,  robust, 
and  efficient  manner. 
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Appendix:  Grid  Size  Adjustment 


Once  the  desired  grid  size  is  determined  based  on  the  truncation  error  as  described  in  sec¬ 
tion  3.2,  it  may  be  adjusted  to  ensure  good  coarsenability.  For  example,  if  the  desired  grid 
is  38  X  28  (grid  intervals  in  x  and  y),  this  can  be  coarsened  only  once  (to  19  x  14)  for  multi¬ 
grid  processing;  since  the  coarsest  grid  is  large,  the  multigrid  solver  is  inefficient.  In  such  a 
situation  it  will  pay  to  increase  the  grid  size  slightly,  e.g.,  to  40  x  28  (which  can  be  coarsened 
twice  down  to  10  x  7)  or  even  to  40  x  32  (which  can  be  coarsened  three  times  down  to  5  x  4). 
We  do  this  as  follows.  Starting  with  the  desired  size  x  Ny  (grid  intervals  in  x  and  y): 

1.  Search  for  a  grid  slightly  larger  than  desired  by  increasing  and/or  Ny  so  they  have 
at  least  one  more  factor  of  two  in  common,  thus  allowing  at  least  one  more  coarse  grid 
for  multigrid  processing.  Accept  the  larger  grid  only  if  it  can  be  strictly  nested  and 
gives  a  predicted  gain  in  efficiency.  Repeat  until  no  gain  in  efficiency  is  produced. 

2.  Test  grid  sizes  slightly  smaller  than  desired  (while  keeping  at  least  as  many  coarse 
grid  levels  as  obtained  in  step  1).  Accept  such  a  grid  only  if  it  results  in  a  substantial 
predicted  gain  in  efficiency. 

To  predict  the  gain  in  efficiency  for  a  proposed  change  of  grid  size,  let  W(;{Nx,Ny)  and 
W^{Nx,  Ny)  represent  the  computational  work  required  to  predict  and  solve  for  tp,  re¬ 
spectively,  for  one  time  step  on  a  computational  grid  of  size  Nx  x  Ny.  We  assume  that 
W^{Nx,Ny)  is  proportional  to  the  number  of  grid  points,  and  thus  to  the  product  NxNy. 
We  also  assume — at  least  for  grids  which  are  sufficiently  coarsenable — that  W^{Nx,  Ny)  « 
u)W(^{Nx,  Ny)  for  some  constant  factor  uj;  empirically  we  find  that  u  ^  2.  Then  the  gain  in 
efficiency  (i.e.,  the  speedup)  in  changing  from  grid  size  Nx  x  Ny  to  Nx  x  Ny  is  predicted  to 
be 

Wc{Nx,Ny)  +  W^{Nx,Ny)  1 

WciNx,  Ny)  +  W^{Nx,  iv,)  ~  a;  -f  1 
To  estimate  we  first  estimate  the  work  required  for  a  multigrid  V-cycle  with  m  grid  levels 
(using  Uc  relaxation  sweeps  on  the  coarsest  grid  and  Uf  on  each  of  the  finer  grids)  as 

Wv{m)  =  Uf(l  +  ^  +  —  +  ---  +  -t-  (19) 

work  units  (one  work  unit  is  the  work  of  one  sweep  on  the  finest  grid),  where  we  have 
neglected  the  work  of  grid  transfers.  The  corresponding  work  for  the  FMG  algorithm  with 
one  V-cycle  per  level  is  then 

Wpim)  =  Wv{m)  -I-  -Wv{m  -  1)  -I - h  ^^_2l^v(2)  -t-  .  (20) 

work  units.  The  number  of  sweeps  on  the  coarsest  grid  is  set  by  the  requirement  that 
(pg)"'  «  (p)""^)  where 


^W4Nx,Ny)  ^  NxNy 
'^WJNx,Ny)  NxNy, 


(18) 
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is  the  spectral  radius  of  the  Gauss-Seidel  iteration  operator  on  the  coarsest  grid  and  p,  =  0.25 
is  the  multigrid  smoothing  factor  for  Gauss-Seidel  relaxation  with  red-black  ordering.  Solving 
for  i>c,  substituting  from  (19)  in  (20)  and  multiplying  by  the  number  of  grid  points  on  the 
finest  grid  yields  the  estimate 


W4N,,Ny) 


N,Ny 

4m-l 


16 ,  _  1  X  4 . 
_(4-i  _  1)  _  _(^ 


l)  +  m 


iog(/^) 

log(pG) 


I 


(22) 


for  the  total  number  of  operations  required  to  solve  for  ^  on  a  grid  of  size  Nx  x  Ny  (here 
normalized  by  the  number  of  operations  required  per  gridpoint  per  sweep).  Combining  (22) 
with  (18)  gives  an  effective  way  to  predict  the  gain  in  efficiency  when  considering  a  new  grid 
size. 
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